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The motion involved in barrier crossing for protein folding are investigated in terms of the chain 
dynamics of the polymer backbone, completing the microscopic description of protein folding pre- 
sented in the previous paper. Local reaction coordinates are identified as collective growth modes 
of the unstable fluctuations about the saddle-points in the free energy surface. The description 
of the chain dynamics incorporates internal friction (independent of the solvent viscosity) arising 
from the elementary isomerizations of the backbone dihedral angles. We find that the folding rate 
depends linearly on the solvent friction for high viscosity, but saturates at low viscosity because of 
internal friction. For A-repressor, the calculated folding rate prefactor, along with the free energy 
barrier from the variational theory, gives a folding rate that agrees well with the experimentally 
determined rate under highly stabilizing conditions, but the theory predicts too large a folding rate 
at the transition midpoint. This discrepancy obtained using a fairly complete quantitative theory 
inspires a new set of questions about chain dynamics, specifically detailed motions in individual 
contact formation. 



I. INTRODUCTION 

A complete description of protein folding kinetics must 
account for the microscopic dynamics of the polypep- 
tide chain. The relaxation of a free polymer is well de- 
scribed by Langevin dynamics with an harmonic poten- 
tial between adjacent monomers and Markovian random 
forces. El This Rouse-Zimm formalism can be extended 
through frictional forces with memory to effectively de- 
scribe chain motions that include many other strictly mi- 
croscopic complexities such as chain stiffness Bel To de- 
scribe the chain motions involved in protein folding, one 
must also attend to the residue specific interactions or 
free energy landscape which guides the protein to its na- 
tive state and governs the probability of different ensem- 
bles of structures described in terms of order parameters. 

For many proteins, the only thermodynamically dis- 
tinct states with significant detectable populations are 
the uafolded ensemble and the native folded conforma- 
tions .□ This observation suggests the existence of a 
free energy barrier separating the stable and metastable 
phases. Like nucleation in phase transitions, the folding 
kinetics for these proteins is largely governed by the prob- 
ability of rare configurations corresponding to the top of 
the barrier, the so-called transition state ensemble. In 
theoretical studies, the transition state ensembles corre- 
spond with saddle-points in a multi-dimensional free en- 
ergy surface as a function of local order parameters. Near 
the transition state, some order parameter fluctuations 
grow to become the folded phase. This is an unstable 
mode of collective motion. Motion in other directions 
tend to stabilize the transition state. The structure of 
the unstable growth mode can be thought of as a locally 
defined reaction coordinate. If there are multiple saddle 



points which are connected, the local reaction coordinate 
for each of these will usually differ. In such cases, espe- 
cially, the local reaction coordinate may not be the best 
way of describing the global transformation of the ensem- 
ble involved in folding. 

Lately, there has been a great effort to elucidate the 
nature of the transition state ensemble and ideal reac- 
tion coptdinate for specific proteins through measure- 
mentjETLfJ theoretical calculations based pa physically in- 
tuitive choices ,EjtLj and simulations Oti3 While many 
aspects of the folding kinetics seem to be understood, 
still some basic issues remain controversial. For example, 
while the structural distribution of the transition state 
ensemble is well predictedjHHl3 the absolute magnitude 
of the barrier height is not well determined. Compari- 
son with experiment is clouded by the uncertainty of the 
prefactor of the folding rate. 

For simple thermally activated transitions, one expects 
the rate to follow the activated law 

k = koe^ AF \ (1) 

where (3 = is the inverse temperature and AF^ 

is the barrier height which contains both entropy and 
energy. For folding, AF^ is expected to be strongly tem- 
perature dependent. The exponential factor is the rel- 
ative probability of finding the system at the transition 
state compared to the metastable minima. The prefactor, 
k , is determined by the microscopic dynamics involved 
in crossing the barrier. For proteins with very rugged 
landscapes, the dynamics, can be described as hopping 
between trapped states£3 In this case the barriers to 
escape traps enter strongly. For proteins with very lit- 
tle frustration that have smoother landscapes, there will 
be little trapping, but the chain motion will be diffusive 
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and depend on chain connectivity, kg depends on the free 
energy surface and the frictional forces such as viscosity 
that determine the stochastic motion of the protein as it 
folds. One estimate proposed for the prefactor relating 
ko to the inverse time it takes to form aijjpical loop given 
by simple polymer theory, ko ~ l/is^O'CJ While quite 
reasonable, this estimate does not involve the shape of 
the free energy surface near the transition state. Thus, 
for example, it is not completely clear how big a loop 
should be considered. 

In Ref. 12] and the preceding paper (henceforth re- 
ferred as [1 77 we presented a variational theory for the 
protein folding free energy surface of a specific protein 
to a known structure. In this approach, the protein is 
modeled as a collapsed, stiff chain with interactions be- 
tween residues in contact in the native structure. The 
reference Hamiltonian is that of a polymer chain inho- 
mogeneously constrained to the native structure. The 
free energy profile for folding obtained from model can 
be described as an activation barrier between the un- 
folded and folded ensembles with fine structure from lo- 
cal minima and saddle-points in the free energy surface 
giving rise to longitudinal ruggedness in the free energy 
profile for a perfectly funneled landscape. In this ap- 
proach, conformational ensembles along the folding route 
are characterized by structural order parameters evalu- 
ated as averages over the distribution from the reference 
Hamiltonian. The focus of the present paper is how the 
microscopic chain dynamics are involved in barrier cross- 
ing within this model. 

Understanding barrier crossing dynamics in protein 
folding requires not only a good characterization of the 
free energy surface and structural ensembles, but also 
a good description of the microscopic motions of the 
polymer backbone. In the Rouse-Zimm model of poly- 
mer dynamics, the motion of the polymer chain is gov- 
erned by the chain connectivity, represented as harmonic 
entropic springs connecting adjacent monomers, and is 
damped, by the interactions of the monomers with the 
solvent .a A more detailed description of polypeptide dy- 
namics acknowledges there are activated transitions be- 
tween gauche andpirans bond angles determined by the 
dihedral potentialE!fLj One way to incorporate the short 
range structural information arising from these barriers 
to internal rotation is through a static chain stiffness that 
defines renormalized harmonic correlations .BE3l3 How- 
ever, a complete description must also incorporate the 
activated nature of these rotational isomerization rates 
in the dynamical correlations. 

Quite early in the study of polymer dynamics, anoma- 
lies found in intrinsic viscosity measurements indicated 
that there was a souEce-af friction independent of the vis- 
cosity of the solvent .lZIlS Kuhn and Kuhn first suggested 
this "internal viscosity" was a consequence of-«iicroscopic 
barrier hopping in the local chain segments £3 This phe- 
nomenological friction found for high polymers may have 
a variety of microscopic origins in addition to internal 
friction due to barrier hopping of many elementary chain 



links. A helpful (though perhaps dated) evaluation of 
different theories can be found found in Ref. |3^. Fix- 
man has shown that simulations for polymers with rigidly 
fixed bond anglcS|-and, bond lengths can be mimicked by 
internal viscosity .E§cil Still other possible sources of the 
measured internal viscosity include the effect, of inter- 
actions with distant monomers in the chain,c3 or other 
inherent noaJinearities neglected in the Gaussian model 
assumption.EJ To our knowledge, the microscopic inter- 
pretation of internal viscosity for artificial high polymers 
has not been definitively settled. It is clear, however, 
that polypeptides will at least exhibit internal viscosity 
from the difficulty of making backbone angle transitions. 
These transitions are roughly 100 times slower than if no 
barriers were present in the Ramachandran potential. 

In this paper, we complete the logical development of 
folding rate theory for minimally frustrated proteins by 
analyzing the microscopic dynamics of chain motions in- 
volved in barrier crossing of the protein folding model 
presented in Ref. [l^ and [I]. We use multi-dimensional, 
non-Markovian Langevin dynamics to describe the bar- 
rier crossing, and identify the unstable modes as the re- 
action coordinates. For the folding rate calculation, we 
use a generalization of Kramer's theoryQ that treats non- 
Markovian dynamics in many dimensions. c3 At the same 
time, the interpretation of a local reaction coepdinate fol- 
lows Langer's theory of the nucleation rates. c£l The de- 
tails of chain dynamics enter the formalism through a 
generalized friction matrix with memory. The friction 
that damps the motion of the chain includes the solvent 
viscosity and an internal friction that depends on the 
length scale (in sequence) of the motion but which we 
take as is independent of the solvent viscosity. The inter- 
nal friction accounts for the activated nature of dihedral 
angle flips in polypeptides. We restrict our attention to 
over-damped dynamics, though it is easy to generalize 
the calculations to include inertial terms as well. 



II. REACTION COORDINATES AND FOLDING 
RATES 

To calculate the barrier crossing rate, one must first 
identify the order parameters that describe the phases 
of the system, and then study their motions. In den- 
sity functional theory, one often considers a density field 
that is slowly varying in space, and expands the density 
functional about a uniform density .lZI This expansion 
yields a free energy consisting of a bulk free energy and 
can be expanded in gradient terms that ultimately rep- 
resent the surface tension between two phases. The den- 
sity then evolves either hydrodynamically or diffusively 
under the influence of the free energy driving force. The 
quasi-equilibrium part of this approach has been applied 
to nucleation of the folded state in proteins starting with 
a mean field free energy functional and defining a global 
native depsity field as a measure of similarity to the na- 
tive stateuS 
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(3F[{ P }]k const + \Y,?¥pJPb 



(8) 



where {r^} specify the configuration of the protein chain, 
and {rf} are the corresponding coordinates of the na- 
tive structure. In the present work, we wish to have a 
description of the folding that is local in sequence and 
hence consider each term in the sum separately. This 
local native density defined as 



P N (r,ri) = 6(r-r;)exp 



2a 2 



a N (r, 



„N\2 



(3) 



where we have relaxed the measure of similarity of the 
native state by replacing S(ri — r N ) by a Gaussian mea- 
sure. 

To connect this density to our folding route calcula- 
tions, it is more convenient to consider scalar order pa- 
rameters than to expand the free energy in terms of field 
variables. Therefore, we integrate p(r, r^) over r, to give 
the scalar measure native similarity 



exp 
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[To simplify notation, we have suppressed the superscript 
N since the native density can still be distinguished 
by the 1-body monomer density, p 1 (v) 1 given in Eq.(I- 
25)]. Following the interpretation of the order parame- 
ters given in [I], we calculate the average native density 
of site i, pi = (p(ri))o, for a given set of variational con- 
straints {Ci} 



Pi [{C}}= Jdvp\{v)p{v) 

= (l + a N G ll )- 3 / 2 exp 



N\21 



3 a N (s 2 -rf) 



2a 2 l + a N d 



(5) 
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that is, pi is a function of {Ci} through the correla- 
tions, Gij, and average position, s^, given by Eq.(I-21) 
and Eq.(I-22), respectively. The free energy can then be 
parameterized by {pi} with -F[{p}] = -F[{C}] where {p} 
is evaluated at {C}. In particular, we denote a saddle- 
point of ^[{C}] by {C*} and the corresponding native 
density by p* = Pi[{C*}], 

We assume Langevin dynamics with memory for the 
native density, {pi{t)}, 



dtPi(t) = - V / dt'p^t-t') 
„■ Jo 



90F[{p(t')}] 
dp] 
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where (3 — 1/ksT is the inverse temperature, and p(t) is 
a generalized mobility matrix related to the random noise 
£i(t) through the correlations (£i(t)£j(t')) = pij(t — t'). 
The form of p,(t) is discussed in the following section, 
here we assume that p(t) is a known function of time. 

To study the dynamics near the saddle-point, we ex- 
pand .F[{/o}] about p* to second order 



where dpi 



p*, and the Hessian matrix 
d 2 PF[{p*)\ 



r* 



dptdpj 



(9) 



has one negative eigenvalue since it is evaluated at a 
saddle-point. In terms of derivatives with respect to the 
variational parameters, f * is the solution of the matrix 
equation 



v- f * dpi dp* _ d 2 pF[{c*}} 



A7 



dGdd 



(10) 



where the derivatives with respect to Ci can calculated 
as described in [I] . 

With this approximation, Eq.(p]) can be written as 



dt(6 Pi (t)) 



(ii) 



Accordingly, the native density correlation functions C(t) 
with components 



Cij(t) = (S Pi (t)S Pj (0)) o , 



satisfy 



d t C(t) 



dt'n(t-t')T*C(t'). 



(12) 



(13) 



In Eq.(|il|), we have subtracted the stationary value 
Pi(t — > oo ) = p*, so that Spi(t — > oo ) = and 
C(t — » oo) = 0. This long time solution is unstable, 
however, because it is a saddle-point of the free energy. 

The unstable mode can be determined by the average 
equation of motionO Denoting the Laplace transform of 
an arbitrary function of time by g(uj) — f^°dte~ u ' t g(t), 
the solution of the Eq.(|l]) averaged over the random 
noise is 



5pi(u) = [wl + £(w)f*].. S Pj (0). 



(14) 



By inverting the Laplace transform, Pi{t) can be ex- 
pressed as a sum of exponentials with time constants 
determined by the poles of [ojl + p(oj)T*^ : the eigen- 
functions p(— n)T* u = nu have time dependence u(t) ~ 
e~ Rt . Assuming that p(uj) is positive definite, there is 
one mode for which k is negative, because T is the saddle 
point Hessian of F[p\: 



u*(t) = u*e 



\k* \t 



(15) 



where 
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A(-«*)F* 



(k* < 0). 



(16) 



Since motion along the unstable mode u* grows exponen- 
tially away from the saddle point, we identify the compo- 
nents of u* as the local reaction coordinate to surmount 
the barrier. Both u* and k* are essentially the unstable 
mode and curvature of F[p], but renormalized by the dy- 
namics of the barrier crossing incorporated in the mobil- 
ity matrix. In Eq.jl^), k* is the multi-dimensional ana- 
logue pfithe Grote-Hynes frequency (in the over-damped 
limit ).e3 

fhe rate for barrier crossing corresponding to Eq.fll] 



isiP 



2tt 



dctr 



det 



1/2 



-/3AF T 



(17) 



where AF is the barrier_ height and f ms and f * are the 
curvature matrices of F[p] evaluated at the metastable 
minimum and the saddle-point, respectively. Eq.(|l7j) 
a genexalizes both the rate-, calculations presented by 
LangeiO and Grote-Hynes it simultaneously accom- 
modates both multi-dimensional diffusion (as in the 
Langer formula) as well as time dependent friction (as 
in the Grote-Hynes formula). 

The ratio of determinants accounts for the entropic dif- 
ferences between the transition state and the metastable 
phase due to fluctuations of the order parameter. These 
fluctuations are like the capillary waves which renormal- 
ize the surface tension in ordinary nucleation. These con- 
tributions would already be included in an exact free en- 
ergy so that the determinants should be absorbed into 
the exponential factor to most .simply keep a consistent 
level of thermodynamic theory £3 Consequently, the ex- 
pression for the rate becomes 



k = e 

2?r 



-/3AF 1 " 



(18) 



where k* is given by Eq. ( |l6| ) . As pointed out in Ref . |5l], 
if the over-counting of entropy were not corrected, the 
ratio of the forward and backwards rate would not equal 
the equilibrium constant as determined by the starting 
free energy functional, k^/kzi ^ e~^^ F2 ~ Fl \ 



III. LINEAR RESPONSE APPROXIMATION OF 
THE GENERALIZED MOBILITY 



assumed to be same as that l^tic} 1 arises from the Mori- 
Zwanzig projected dynamics.E3E3 We begin by outlining 
this formalism as it was applied in Zwanzig's classic paper 
on generalized Rouse dynamics for a polymer described 
by an arbitrary potential.0 A closed form of the mobil- 
ity matrix is then obtained by approximating the local 
dynamics with those of the constrained polymer Hamil- 
tonian, H$. 

Consider a polymer described by the potential U(R) 
with monomer positions R = (r^, . . . , r„). The probabil- 
ity density of monomer positions, ^(R, t) is assumed to 
evolve according to the Smoluchowski equation 



<9 t *(R,<) = 2?*(R,t) 



P = ^V l -A, e - /3C/(R) V,e^ R ) 



(19) 



(20) 



where V.£ = d/dr^ and D is the diffusion matrix. In this 
section, we denote the equilibrium averages by (...) = 
JdR . . . * eg (R) where 



* eq (R) = e -^ R > 



(21) 



Now consider a set of dynamical variables that are 
functions of the monomer positions, {Ai(R)}. These 
variables will be chosen to be the order parameters {p} 
to describe the barrier crossing dynamics discussed in 
the previous section, but for now {^4} is arbitrary. The 
equilibrium time dependent correlation functions of 



SAi =Ai- (A,) 



are denoted by 



cfJt) ee (6Mt)SM0)). 



(22) 



(23) 



The time dependent correlations can be formally ex- 
pressed as 



<#(<) 



dRSA l e m (6A j ^ eci ) 



dR^^SAie^SAj. 



(24) 



In the second line of Eq(24|), we have introduced the ad- 
joint operator C defined by its action on an arbitrary 
function B 



The time evolution of the native densities {Spi(t)} de- 
pends on microscopic dynamics of the monomer positions 
{ri(t)}. The chain dynamics, in turn, are incorporated 
into the effective mobility matrix p(t) in the relaxation 
equation for 5 pit) (Eq.flll])). To determine p,(t) consis- 
tent with the definition of the native densities this theory, 
we consider the chain dynamics near the saddle point to 
be governed by the Hq with constraints {C*}. The re- 
sulting mobility matrix for native densities in Eq.(|ll|) is 



WS> cq B = V cq £B, (25) 

or explicitly, 

C = £ e^WViDije-WW ■ V,, (26) 

ij 

The meaning of the operator notation in Eq.(^4j) is 
defined through the propagator (or Green's function) 
P(Xi;Y): 
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<#(t) = JdXdYSA t (X)SA 3 (Y)P(Xt;Y)^ cq (Y), (27) 

where (formally) 

P(X <; Y) = e I3 '(5(X - Y), (28) 

i.e., P(X t;Y) satisfies the Smoluchowski equation, 
d t P = VP, with initial conditions P(Xt = 0;Y) = 
S(X - Y). 

An alternative equation can be derived for the corre- 
lations through the projected dynamics of Eq.([l9|). Fol- 
lowing Ref. ^], we define the projection operator through 
its action on a arbitrary function B 

VB = J2 5A i r tj( 5A i B ) ( 29 ) 

ij 

where (r A ) _1 is the matrix of static correlations 

(r A )- x = (SAiSAj). (30) 

Using standard projection operator techniques^ it 
can be shown that-the correlations obey the generalized 
Langevin equation!] 

d t c A (t) = -n A r A c A + f dt'K A {t - t')r A c A {t'), (31) 
Jo 



where 



-(AiCAj) 



kl 



[WkAi-DkrWiAj) (32) 



and K (t) is the memory kernel 

K A (t) = (Ai£(l - P)exp[(l - V)Lt]{l - V)CA 3 ). (33) 
Eq. ( |3l| ) can be written in the form of Eq. ( |l2| ) 



d t C A (t) = - / dt'p, A {t - t')r A C A (t'), (34) 
Jo 



where 



fi A (t) = 2n A S(t) - K A {t) 



(35) 



is the effective time-dependent mobility matrix. As sug- 
gested by this notation, the static correlations (r A ) _1 
can be thought of as arising from an effective harmonic 
free energy j3F = const + (1/2) £\. JAT^Aj . 

Eq.(^) is an exact result, equivalent to Eq.(|27j). The 
linear response formulation avoids explicit use of the 
Green's function which is generally unknown for an ar- 
bitrary potential. On the other hand, it is very difficult 
to calculate the memory kernel explicitly from Eq.(^) 
since the projection operators are very awkward to ma- 
nipulate. 

If we choose {A} to be the monomer positions, then 
Eq.(|3l|) becomes an equation for the monomer correla- 
tion. If one neglects the memory term, the monomer oa: 
relations are said to obey optimized Rouse dynamics 3 



This is equivalent to having an effective pre-averaged dif- 
fusion matrix, f2 = (D), and a Gaussian approximation 
to the chain through the static correlations. [Bixon and 
Zwanzig first presented the stiff chain connectivity matrix 
[Eq.(I-4)] in a paper that employs the optimized Rouse 
dynamics formalism]. One usually must resort to ap- 
proximate methods to include the memory kernel in the 
analysis. A clear presentation of these methods applied 
to one-dimensional problems is given in Ref. ^4]. 

One way to study the memory function for monomer 
correlations is through a truncated Mori continued frac- 
tion.E3 This becomes rapidly more cumbersome as the 
order of the expansion is increased. A more convenient 
expansion of the memory kernel is to first expand the 
Green's function in Eq.(|2^) in the eigenfunctions of C, 
and then expand these eigenfunction in a finite basis 
(similar _tp_what is done in quantum chemistry calcu- 
lations). EjI£3 The basis functions must be chosen with 
some care for accurate results. Furthermore, in both the 
continued fraction and eigenfunction expansion methods 
require equilibrium averages that are difficult to com- 
pute for an arbitrary potential, requiring simulations in 
general. Nevertheless, this is a useful method since the 
simulation time required for the equilibrium averages can 
be much shorter than would be required to simulate the 
time-dependent correlation functions themselves .E3 

To connect to the barrier crossing calculation for pro- 
tein folding, we need to study the dynamics of the native 
density, {A} — {p}. For a general potential, calculating 
the correlation functions (5pi(t)5pj(0)) involves the same 
technical difficulties that arises for the calculation of the 
monomer correlations. We avoid these complications by 
approximating the potential by the constrained polymer 
Hamiltonian, Ho. Not only is this is a reasonable choice 
given the interpretation of the order parameters as aver- 
ages over Hn, it also allows us to calculate the correlation 
from Eq.(j27|) directly. 

For an harmonic potential, the propagator in Eq.(|27]) 
can be solved analytically^ The Gaussian form of p(r^) 
allows the integrals in Eq.(p7|) to be performed to give 
the native density correlations (see Appendix for details) 



(pi(t) Pj (0)}o = (detM^i))" 3 / 2 
-3a N 



x exp 



2a 2 



jf r Mutt)- 1 Jt 



where 

Afy(t) = 

and, Jij is the 2-component difference vector 



a N Gy(i) (l+oFGy) 



■i 



(36) 



(37) 



(38) 



In this expression, we have introduced the time correla- 
tion function of monomer positions 
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G ij {t) = (6r i (t)-8r j (0)) /a 2 



(39) 



IV. CHAIN DYNAMICS 



The chain dynamics are contained in Gij (t) , and the pre- 
cise form depends on the diffusion matrix Dij in Eq.(20). 
We postpone further discussion of Gij(t) until the next 
section. 

The correlations in Eq.(p7|) for {A} = {p} are written 

as 



Cy(t) = (5pi(t)6p(p)) = ( P i(t) Pj (0)) - PiPj (40) 



where the first term is given by Eq. (J36|-p8[) and the second 
term is found from Eq.(|^). 

Since approximating the chain potential with Hq per- 
mits an exact solution for the correlation functions, it 
may seem as if the discussion of the projected dynamics 
was an unnecessary detour. However, what we are really 
trying to calculate is the effective mobility matrix needed 
to evaluate the folding rate prefactor in Eq.([l6|). 

From Eq.(34), the native density correlations obey the 
projected relaxation equation 



d t C{t) = - / dt'p p {t - t')T p C(t'), (41) 
Jo 

where T p are the inverse of the static correlations 



(dpiSpj)o = Cy(0). 



(42) 



Eq. (p5[) with {A} = {p} gives p p (t), in terms of pro- 
jection operators, but we can also deduce an alternative 
expression for the mobility since C(t) is known. Solving 
Eq.([4l"l) for p p (t) by Laplace transforms gives 



A P H = c(o)dH _1 c(o) - wC(o) 



(43) 



To close this expression for the mobility, we only need 
the Laplace transform of C(t) which can be calculated 
numerically. 

As discussed previously, Eq.(f4l|) is consistent with a 
Gaussian approximation to the free energy 



(3F P = const 



-Yt p 



SpiSpj 



(44) 



analogous to the saddle-point expansion of given 
in Eq.(||). However, even though T p can be evaluated at 
the saddle-point constraints in {C*} in Hq, F p is stable 
against all fluctuations away from dpi = 0. This is in 
contrast to the curvature matrix T* that has one unstable 
direction. 

Nevertheless, by approximating the mobility in Eq.(|l2]) 
by u(t) ~ fJ- p {t) we can solve the eigenvalue equation 
Eq. (|l6|) to identify the unstable mode u* and calculate 
the prefactor k*. The last remaining unspecified quan- 
tity for the theory is the time dependent monomer cor- 
relations G(t) which is the subject of the next section. 



In this section we specify our model for the chain 
dynamics, and calculate the monomer pair correlation, 
G(t), with the approximation to the chain potential 
U(R) = Hq. This Gaussian approximation simplifies the 
true microscopic potential that governs conformational 
changes of the protein chain, and by itself ignores ex- 
plicit interactions. We attempt to capture some effects 
of the microscopic dynamics by including frictional forces 
that are independent of the solvent viscosity, i.e., internal 
viscosity. 



A. Monomer Pair Correlations 

Starting from the definition of pair correlations for the 
Smoluchowski equation [Eq.(p7|)], it is easy to derive an 
equation for the relaxation of the monomer correlations 

dtiSr^t) • <5r,(0)) = <V fc tf({r(t)}) . 5r .( )), (45) 

k 

where the friction coefficient is related to the diffusion 
coefficient by 



(46) 



In Eq.(45), we have assumed that 7 is independent of 



{r}, but may still depend on the sequence index. The 
model for the friction matrix is specified shortly, for now 
we define a dimensionless friction matrix, 7: 



7 = 7o7 



70 = 67ra off 7/, 



(47) 



where we have used Stokes law to define 70 with the 
solvent viscosity 77 and typical monomer van der Waals 
radius a e g. 

For a given point along the folding route (e.g., a transi- 
tion state) , we assume that the chain dynamics can be de- 
scribed by the reference Hamiltonian. [Eq.(I-13)]. Physi- 
cally, Hq corresponds to a polymer confined to the native 
structure by harmonic constraints with spring constants 
{3/2a 2 Ci}. These constraints act as a non- uniform ex- 
ternal field, controlling the fluctuations about the native 
structure. Taking the gradient of Hq, and using the def- 
inition [Eq.(I-15)] of the average monomer position, Sj, 
gives 



dH _ 3k B T 
dvi a 2 

3k B T 



(ch). 



(o) Sr . 



(48) 



where is given in Eq.(I-17) and Sr^ = — Sj as usual. 
Substituting this into Eq.(f45|) gives the matrix equation 
for the correlations Gij(t) = (6ri(t) ■ Srj(0))/a 2 , 



G 



where 



d t G(t) = aT^Git) 
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(50) 



sets the time scale for the relaxation. 

Since Eq.(|49"|) is linear, it can be solved by transforming 
to normal modes by a similarity transform: 



Q 



^ lr r 1 r (0) ]Q = diag{A p } ! 



(51) 



where the columns of Q are the right eigenvectors of 
py -1 IH°)]. Even though 7 and r^ ^ are symmetric, their 
product is not symmetric in general. Nevertheless, Q 
can diagonalize 7 and I^ ' separately (though not by a 
similarity— .transformation since Q is not orthogonal in 
general) S3 With 



Q T iQ 



diag{f p } 



(0) 



we have the relationship 



Ar, 



r] p /v p . 



Q = diag{r, p } (52) 



(53) 



The correlations between monomers, expanded onto the 
normal modes is 



G ij(t) =^2Q P iQjp/v P exp(-aX p t). 



(54) 



We note that equal time correlation functions are the 
equilibrium correlations Gij(t = 0) = Gy = [r^ )] -1 . 

The reference Hamiltonian has two features that dis- 
tinguish this chain model from the nearest neighbor 
Rouse chain: chain stiffness, parameterized by g, and 
constraints to the native structure, specified by {C}. Wc 
consider the effect on the pair correlations of these two 
contributions separately, to give a point of reference for 
the model. To simplify this illustration of the dynam- 
ics, we consider a (free-draining) diagonal friction ma- 
trix, 7y = Sij. Then, v p = 1 which implies that the 
relaxation rates are X p — rj p and the normal modes are 
the eigenvectors of the chain connectivity r( ch ). 

For an unconstrained polymer (Cj = 0), the dynamics 
of the stiff, jfceely- rotating chain has been studied quiet 
extensively BE^Te 2 ] Referring to the definition of T [Eq.(I- 
4)], Q differs from that for the well known rouse modes 
Q R 



[Q*fK R Q R = diag{<} 



(55) 



because of the boundary term, g2-I(l — g)A, which is only 
important for very stiff chainsHj Thus, unless g ~ 1, 
Q k Q R for the free-draining, unconstrained polymer: 



cos 



n 



Ti/ p = 4 sin 2 



np 
2n 



(56) 



where the mode index takes values p = (0, . . . , n 
The relaxation spectrum is then approximately 



Ap — Vp 



9 R 
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B 



(57) 



where the final term is due to the confinement term in 
Eq.(I-9). The second term of Eq. (|57|) causes the short 
wavelength modes (high mode index p) to relax faster as 
the chain stiffness increases. In contrast, the relaxation 
rates for the long wavelength modes (low mode index p) 
rates are suppressed by increasing chain stiffness. The 
number of modes with slowly relaxing rates that are sup- 
pressed diminishes as the chain stiffness increases; for the 
limiting stiffness corresponding to a rigid rod, g 1 , the 
relaxation rates for all modes except one diverge. 

When the polymer is constrained by a non-uniform 
pattern of nonzero constraints {Ci}, the normal modes 
are not the Rouse modes of a free chain. Fig|l] shows 
the four slowest modes for a chain with stiffness param- 
eter g = 0.8 (persistence length I = 5 a) and various 
constraints. In addition to the modes for the uncon- 
strained chain shown in Fig.[l]a, the modes correspond- 
ing to selected stationary points along the folding route 
of A-repressor are also shown in Figjljj-f. The slow- 
est modes for the denatured state (Fig.[I]b) differ only 
slightly from those for the unconstrained chain (Figjlja) 
given by Eq.(|56]). Referring to the fluctuations shown 
in Fig.I-9b,the structure at the transition state TS± can 
be roughly described by saying that residues in helices 
H4-H5 (indices 53-80) are localized, while the residues 
in helices H2-H3 (indices 27-47) have the largest fluctu- 
ationsEJ . These constraints are reflected in the structure 
of the slowest relaxation modes show in Fig.^jc. The mo- 
tion of the more localized residues is shifted to modes 
with higher relaxation rates so that these residues are 
absent in the slowest relaxation modes. In addition, the 
mode with the slowest relaxation (solid line) is no longer 
a uniform translation mode, but involves motion of the 
most unconstrained resides in helices H2-H3 (indices 27- 
47). The other modes are much the same as those of 
an unconstrained but shorter chain with one end fixed 
at the end of helix H4. Similar correspondence between 
the mode structure and the localization of residues in the 
transition state ensemble can be seen for the constraints 
corresponding to TS3 and TS4 in Fig.|l|d-e where the 
modes become increasingly more localized in sequence 
due to non-sequential constrained residues. The normal 
modes for the native state (Fig.[j]f) are more localized 
with the slowest relaxation involving localized regions 
with the largest temperature factors. 

The mode spectra (relaxation rates) for the free chain, 
constraints corresponding to transition states, and the 
native state are shown in FigJ^. As can be seen in this fig- 
ure, while the spectrum for the short wavelength modes 
(large mode index) are relatively unaffected, the relax- 
ation rates for all the long wavelength modes (small mode 
index) become faster as the protein takes on more struc- 
ture, reflecting the relatively rapid motion toward the 
native structure as the constraints provided by the na- 
tive structure become stronger. 
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B. Internal friction 

One can expect that Eq. ( f49| ) with harmonic forces can 
capture the correct physics of the low frequency spec- 
trum of the chain dynamics, but, as in the Debye the- 
ory of crystals, the description will not be correct for 
small length scales which depend on microscopic details. 
For many polymer physics applications, even though the 
polymer has a broad range of relaxation rates, the ex- 
perimentally observed dynamics is dominated by the the 
long wavelengths (in sequence) motion corresponding to 
the slowest modes, de Gennes calls the independence 
of the long wavelength modes on the micccecopic con- 
formational changes the "Kuhn Theorem" e3 Why then 
is it interesting to incorporate internal viscosity (arising 
from local activated isomerization) into a model of bar- 
rier crossing dynamics? Because we do not know a priori 
the scale of the motions involved in locally crossing the 
transition region. 

To be more specific, we consider internal friction for 
the following two reasons: (1) We have seen that the 
nonuniform constraints to the native structure affects 
the structure of the normal modes moving them to in- 
volve more local motions. It is unclear then how the 
Kuhn theorem (which is a statement about the delo- 
calized continuum limit) applies here. (2) We are ulti- 
mately interested in studying the barrier crossing rate 
for protein folding, through the Laplace transform of the 
native density correlations [Eq.(|3^)] which mix different 
Rouse modes. A correlation function similar to the the 
native density correlation function arises (for an uncon- 
strained polymer model) in the Wilcmski-Fixma^thc- 



ory of diffusion-limited chain closure reactions.E§E3 In 
analyzing this theory, Doi showed that for end-to-end 
monomers, the Laplace transform evaluated at small fre- 
quencies is actually dominated by the integral at short 
times, i.e., all modes contribute to the iategral, not just 
those with the longest relaxation times. E3 

The friction matrix in Eq.(]49|) 7 = 707 is written a 
sum of two contributions 



As previously discussed, some internal friction is due 
to the neglect of the activated character of the individ- 
ual conformational transitions of the backbone angles. 
While there are other ways to incorporate this time scale 
into a model of chain dynamics perhaps the most sim- 
ple correction to the Gaussian description is through the 
friction matrix. In accordance with the Kuhn theorem, 
the long wavelength motion should have the usual sol- 
vent viscosity dependence being largely unaffected by the 
microscopic barrier crossing; on the other hand, short 
wavelengths should tend to relax more slowly, due to the 
"dynamical stiffness" caused by the slow activation rates 
for dihedral angle changes. Since the dihedral angle po- 
tentials are quite steep these will have a much weaker 
viscosity dependence. 

In the model presented in Ref. |7(], Bazua and Williams 
derived an effective frictional force on the i th bead which 
acts on the relative velocities of adjacent beads: 



(60) 



where the superscript denotes the projection of the forces 
along the bond directions. A common approximation is 
to write the projected velocities as a sum of the difference 
between the full velocity vectors {r} and. the difference 
between angular velocity of the beads.Eil The angular 
component is important for inherently anisotropic phe- 
nomena such as flow birefringence, but for our purposes 
we can neglect the angular rotatkm^omponent and write 
the internal friction in the form&ErEj 



(61) 



where K R denotes the Rouse connectivity matrix [Eq. (I- 
5)], and we have used dt(Sri) = dtTi since the average 
monomer position is constant for fixed constraints. 

The total friction matrix is the sum of the (free- 
draining) diagonal solvent friction and the internal fric- 
tion 



Hi 
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(58) 



7o (Sij 



t-7,/70^) 



(62) 



where the first term is proportional due to the solvent 
viscosity, and the second term accounts for the "internal 
viscosity" and is taken to be independent of the solvent 
friction. A typical choice for the effective solvent friction, 
is the Oseen tensor from the hydrodynamic mobility .El 
This non-local friction mediated by the solvent damps 
the motion of one bead due to the motion of another 
distant beads. It is not difficult to use standard preaver- 
aging approximations to incorporate these hydrodynamic 
interactions. For simplicity, we neglect this aspect of the 
problem and consider the free-draining approximation 



7 S - 



where 70 is defined in Eq.(|47 



(59) 



which defines the dimensionless friction matrix 7 in 
Eq.©. 

The spectrum of relaxation rates with nonzero inter- 
nal friction can be illustrated by considering an uncon- 
strained polymer chain (Ci — 0). From Eq.(p3|) the mode 
relaxation rates for this case is approximately 



(63) 



where -q are the eigenvalues of T^ ) (Eq.(|57|))and rjp are 
the eigenvalues of the Rouse matrix K R . This equation 
shows that the internal friction leaves the rates of the 
long wavelength modes (small rj R ) relatively unchanged; 
this is the Kuhn theorem. The short wavelength modes 
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(larger rj ), are suppressed by a non- vanishing internal 
friction. The crossover between the two behaviors, hy- 
drodynamic solvent limited versus conformational inter- 
conversion limited, is determined by the value of 77/70- 
For constrained chains, the effect of internal friction on 
the relaxation modes is different. The internal friction 
damps rapidly oscillating features of the modes which 
even occur in the relaxation modes with long relaxation 
times since they are not identical to Rouse modes. (Note 
the modes illustrated in Fig.0, which are also different 
from Rouse mode were calculated for 7/ = 0). Neverthe- 
less, for the constraints considered here, this seems to be 
a minor effect, and the relaxation modes still follow the 
general behavior of Eq.(|63|). Fig.|| shows the relaxation 
spectrum with zero and finite internal friction for the un- 
constrained chain and for the constraints corresponding 
to the transition state TS2 of the variational calculation 
for A-repressor. The spectrum corresponding to the other 
transition states is similar: the relaxation rates for the 
slow modes are relatively unaffected, whereas the faster 
modes have suppressed relaxation rates. 



V. FOLDING PREFACTOR EXAMPLE: 
A-REPRESSOR 

The average folding route for the fast folding vari- 
ant of the A-repressor protein presented in [I] is char- 
acterized by sets of constraints {C} corresponding to lo- 
cal minima and saddle-points of the variational free en- 
ergy. For each saddle-point, the pair correlations G{t) 
are calculated through the normal mode expansion given 
by Eq.(|54|) with the friction matrix given by Eq.(62). 
These monomer correlations determine the native den- 
sity correlations C(t) through Eq.(^0|) [with Eq.(||) and 
Eq.(|36|)1, and its Laplace transform C(us) is taken numer- 
ically. With the mobility matrix p,(uj) given by Eq.(|43|) 
and the curvature matrix of the free energy T* given by 
Eq. (|l0|) , the eigenvalue equation Eq.(|l6|) is then solved 
numerically giving the reaction coordinate and folding 
rate prefactor for that transition state. 

In this example, we consider the folding routes obtain 
for a chain with persistence length of polyalanine, I — 5 a 
(g = 0.8). The native density {p} also depends on the 
parameter a N controlling the width of the Gaussian mea- 
sure to the native structure. This width should be small 
enough to characterize the native structure accurately, 
but the appropriate value is also limited by the resolu- 
tion of the model (approximately equal to the root mean 
square (RMS) fluctuations of the bond length a). For 
large a N (small width), {p} has low values even for the 
native state since the average positions are not precisely 
equal to the native structure at finite temperature; for 
cv N = 1 the average value of pi evaluated at the native 
state is 0.3. For the native densities shown in Fig.[|with 



for a N . For smaller a N (larger widths), the separation 
between {p} evaluated at the native and globule state is 
more distinct, though unfolded regions of the transition 
states have native densities closer to the native values. 
The dependence of the prefactor on the width a N is con- 
sidered below. 

The local reaction coordinates (unstable modes) for the 
four transition states along the folding route are shown in 
Fig.^. These unstable modes are quite understandable in 
terms of the structure inferred by the temperature factors 
in Fig.I-9b.The unstable mode of TS\ (Fig.^a) is rather 
delocalized, consisting of residues in helices H4-H5 and 
modulates with a wavelength approximately equal to the 
persistence length l/a = 5; the partial localization of he- 
lix HI indicated by the temperature factors is not a part 
of the normal mode. The mode for TS2 (Fig||b) has a 
significant component of formation of helices H4-H5 as 
in TS\, but the localized region around a residue in HI 
is a much larger component. The unstable mode at TS3 
is also localized, but consists mainly of a few residues 
separated in sequence in helix HI and near helix H4. Fi- 
nally, the unstable mode at TS4 is delocalized consisting 
of residues in helices H2-H3 which are the last helices to 
form. The features of the reaction coordinates resemble 
the unstable modes of the free energy (f), though the 
relative amplitudes of the peaks in a given reaction co- 
ordinate are altered by the chain dynamics through the 
mobility p(to). 

Since the local reaction coordinates along the average 
folding route are very different, one from each other, 
a reduced description of the kinetics based on a single 
reaction coordinate can not be used when dynamics is 
examined at this level of detail. Attempts to find per- 
fect global dynamical coordinates in lattice simulations 
of perfectly tunneled models are consistent with thisH3'E3 
In the present description, each dynamical reaction co- 
ordinate applies in the vicinity of one saddle-point. The 
changes of the local reaction coordinate along the folding 
route are intimately related to the fine structure of the 
free energy profile. Notice that the effects on predicted 
reaction rates are minor, however. Surmounting the high- 
est barrier gives the rate within a factor of two (vide 
supra). Calculations of one-dimensional barrier crossing 
with shallow local minima (intermediates) also show why 
having a perfectly good local detailed description of the 
free energy profile need not yield a strictly defined global 
reaction coordinate.^ At the same time, a much cruder 
reaction coordinate can function perfectly well for pre- 
dicting rates if the_appropriate configurational diffusion 
coefficient is used.O 

The prefactor ko — |k*|/27t for each transition state 
is plotted as a function of the internal friction 71/70 in 
Fig|| When the unstable curvature of the free energy 
is large, there is less re-crossing at the top of the barrier 
leading to a larger prefactor. For small 77/70 the prefac- 
tors increase in the same order as the corresponding cur- 



a = 0.5, the average of pi at the native state increases to 

Pi is 0.5. We consider this to be a reasonable value to use vatures of the free energy (r ): Jtf> < k£> < k < k, 
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The relative values of the prefactors at different tran- 
sition states depend on the chain dynamics through the 
mobility matrix p,(oj), however, and do not follow the pro- 
portions of the corresponding free energy curvatures. For 
small values of 77/70, the prefactor is proportional to the 
bare diffusion coefficient Dq = ksT/^Q. The prefactors 
decrease with increasing 7/ / 70 because of the suppressed 
relaxation rates of the fastest modes. This turnover oc- 
curs around 77/70 = 0(1), though it is different for each 
transition state. 

To put these results in laboratory units for polypep- 
tides, we need to specify the monomer diffusion coef- 
ficient Dq and the internal viscosity parameter 77/70- 
We estimate the typical hydrodynamic radius in Eq. fl47| ) 
to be the sum of the monomer spacing a = 3.8A and 
the average van der Waals radius of a typical side chain 
Q-vdw ~ 3 AEj et ff w 7A. Using the viscosity of wa- 
ter at room temperature r/o = 1 centipoise gives the a 
monomer diffusion coefficient D^ 2 ° = 3 x 10~ 6 cm 2 /s. 
We estimate the internal viscosity by the elementary 
bond flipping rate of the polymer backbone. Calcula- 
tions of the (3 — > a conformations of di-alanine give a 
value for r^t ~ 10ns _1 . E3 One can expect this rate to 
be reduced when the flipping transition occurs in a long 
polypeptide chain since the neighboring monomers must 
adjust to some extent in ordes-tp minimize the motion of 
the chain through the solvent E3 Judging from studies on 
model chains one expects about a factor of three slower, 
we estimate a typical value to be citejs:eh:80,eh:js:82 
= 3 ns -1 . It would be nice for simulators to pin 
this number down for the peptide isomerization rate in 
long chains. It is obvious that the relaxation rates de- 
scribing the chain dynamics should not exceed this ele- 
mentary rate. Associating the internal friction with this 
isomerization rate, we set 77/70 so that the free chain 
relaxation rates of the fastest mode equals r^ 1 . For a 
chain with persistence length I = 5a and monomer diffu- 
sion coefficient D^ 2 ° , we find 77/70 = 9.1. 

In Fig.0, the prefactor fco is plotted against the inverse 
monomer diffusion coefficient Dq 1 in laboratory units. 
The internal friction is set to the fixed value 77 = 9.I70 
for 7o/fcsTf = (D^ 20 ) -1 . The prefactor is proportional 
to the solvent viscosity (D^ 1 ~ 770) for large viscosities 
and saturates at low viscosities. This turnover is remi- 
niscent of the viscasity dependence of the prefactor from 
Kramer's theory JU but here the turnover is driven indi- 
rectly by the internal friction of the chain rather than ex- 
plicit inertial terms reflecting an energy limited regime. 
The prefactor in pure water, indicated by the dashed 
line in FigJ^, is still in the linear regime, but near the 
turnover. The prefactors of the different transition states 
for Do = D^ 2 ° (given in Table are on the order of 
10/7S -1 . As a point of reference, the unstable mode fre- 
quency for TS2, \k*\ = 6.9 x 10 -3 <t, is approximately 
equal to the relaxation rate of the p = 6 normal mode of 
a free chain of n = 80 monomers and persistence length 
I = 5 a (A p=5 = 7.2 x 10~ 3 ). The normal mode for this 
motion has a node to node period of 20 monomers, twice 



the length used-in the estimates based on the smallest 
loop closure. oEj 

The folding rate includes the activation barrier (cal- 
culated with this model in [I]) as well as the prefac- 
tor. The calculated thermodynamic and kinetic quan- 
tities needed for the folding rate are given in Table || for 
two thermodynamic conditions. Adjusting the energy 
scale of the contacts eo/fcsT fixes the equilibrium con- 
stant K = kf/k u (the ratio of the folding and unfolding 
rates) , so that In K gives the difference in free energy be- 
tween the globule and native state (in units of fcgT). For 
the very stabilizing condition \nK — 7.8, the prefactor is 
fco = 2.4/is -1 and the barrier is AF' jksT — 3.4, giving 
a folding time of 1/k = 13/iS. This agrees well with the 
experimental rate 1/k = i3/is at this stability (extrapo- 
lated to 0M denaturant)JZ3 At the transition midpoint 
InK = 0, the calculated prefactor is 7.6/is, and with the 
barrier height AF^/fe^T = 5.1, we get a folding time of 
1/k = 23/7S which is considerably faster than the mea- 
sured folding time of 1/k w 10ms. O This suggest the 
microscopic calculation for the perfectly funneled surface 
under estimates the barrier by about 5fc^T at the folding 
transition temperature. 

The calculated inverse prefactor for A-repressor at the 
transition midpoint 1/fco ~ 0.1/is and under highly sta- 
bilizing condition 1/fco ~ 0.5/zs are botivJaster than the 
estimated fastest possible folding timeE-3c3 r w 1/is. The 
latter is the time it takes to form a typical contact as 
determined by the fluorescence or chemical quenching 
rate of a pair of residues separated by a typical loop 
length. It is very difficult to imagine that folding can 
occur on times faster than contact formation, and the 
introduction of this speed limit is important conceptu- 
ally as a reference timescale for protein folding dynam- 
ics. It is tempting to identify the speed limit with the 
folding prefactor since the fastest rate is obtained from 
Eq.(Q) for a vanishing free energy barrier. But this is 
not precisely correct since the separation of time scales 
assumed in activated kinetics is not true for small barri- 
ers. As seen in Table 0, reducing the activation barrier 
(by stabilizing the native state, for example) changes the 
prefactor since the saddle-point curvature of the free en- 
ergy changes. If the prefactor at both stabilities in this 
example were 1/fco = 0.1/iS, the folding time with the 
reduced activation barrier would be a factor of approxi- 
mately 6 times smaller than at the transition midpoint, or 
1/fc » 4/xs, which approaches the speed limit. The pref- 
actor is found to be approximately three times smaller 
than at the midpoint, however, resulting in a larger fold- 
ing time 1/fc ss 12 fj,s. Pushing the conditions to the limit 
of a vanishing barrier, it could happen that the barrier is 
very low but still has a large unstable mode free energy 
curvature, resulting in a calculated folding time some- 
what fasteij-tlia|n the speed limit. For such near down- 
hill folding JI§Ej however, the folding would not be single 
exponential process dominated by barrier crossing, but 
involves kinetics more directly connected to the time for 
contact formation than the time reflected in the apparent 
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prefactor. 

It is interesting to consider the dependence of the pref- 
actor on the chain dynamical parameters since the values 
for the monomer diffusion coefficient and the isomeriza- 
tion rate in a polymer chain are not known precisely (see 
next section). Conservatively, the values we used are 
probably accurate to within a factor of two or three. As 
long as the internal friction is not too large, the monomer 
diffusion coefficient is the more important parameter for 
the absolute prefactor (in s^ 1 ); since the time scale for 
the monomer relaxation is 3D /a 2 , the absolute prefac- 
tor is proportional to this scaling. Table [II gives the 
prefactor from TS2 using the free alanine diffusion coef- 
ficient and the di-alaninc isomerization rate. As shown 
in Table |l|, the prefactor varies by about a factor of 
three when D^ 2 ° is the diffusion coefficient of a free ala- 
nine molecule. This is the same factor as the ratio of the 
monomer diffusion coefficient, though it also depends on 
the value of the internal friction (i.e., t^*) as well. 

The value of the Gaussian width a N defining the na- 
tive density is a parameter of the model, not a measurable 
quantity, so we must also consider the sensitivity of our 
results to the value of a N . As shown in Fig.||, the pref- 
actors fco are approximately linear in a N , increasing by a 
factor of about two for TS3 and less than a factor of unity 
for the other transition states as a N varies from 0.25 to 
1.0. The different slopes are reflect the different increases 
in curvatures of the free energies as the Gaussian width 
narrows. This degree of uncertainty seems reasonable, 
but the different slopes suggests that the dependence is 
a delicate issue, depending on the detailed structure of 
each transition state. 



VI. CONCLUDING REMARKS 

In this paper, we derived a microscopic theory for the 
dynamics of protein folding barrier crossing when non- 
native contacts and trapping effects can be neglected. 
The variational theory presented in [I] characterizes the 
average folding routes as connected minima and tran- 
sition states. These routes are specified through varia- 
tional parameters constraining the protein chain to in- 
homogeneously order about the native structure. The 
quadratic reference Hamiltonian from this theory along 
with forces from solvent and internal friction govern 
the chain dynamics near the saddle-points. The barrier 
crossing dynamics is calculated from a generalization of 
both Grote-Hynes and Langer rate theories to include a 
frequency dependent friction in multi-dimensions. The 
memory friction itself is approximated by linear response 
theory and exploits the harmonic potential of the model. 
The theory gives local reaction coordinates and folding 
rate prefactors for specific proteins with known native 
structure. 

For the calculated average folding route of A-repressor, 
we find a folding prefactor fco ~ 10 /is -1 at the transition 



midpoint and fco ~ 1/^s -1 under highly stabilizing con- 
ditions. Using the free energy barriers from this model, 
the calculated folding rate agrees well with the measured 
rate extrapolated to zero denaturant, but not at the tran- 
sition midpoint. Considering other reasonable values for 
the parameters in the model suggests that the calculated 
prefactor may be off by a factor of about three. These cal- 
culations also suggest that the folding rate prefactor near 
water solvent conditions scales linearly with solvent vis- 
cosity. This dependence agreesj&iith experimental studies 
on other two state proteinsEjEa and simulations of off- 
lattice model proteins£3 Nevertheless, these conditions 
are near a turnover regime for which the internal friction 
of the polymer chain dominates, ultimately causing the 
prefactor to saturate upon decreasing solvent viscosity. 
Thus folding of polymers with other backbones may be 
quite different. 

The discrepancy between the calculated and measured 
folding rate at the transition midpoint (Tf ) suggests that 
the model underestimates the activation barrier or over- 
estimates prefactor for this stability. Since the model has 
many simplifications that could influence both the ther- 
modynamics and the dynamics, it is impossible without 
further investigation to identify which primarily needs 
improvement. It may be combination of both. 

One aspect of the present model is that the native state 
is quite a bit more flexible than observed in the labora- 
tory, with RMS fluctuations on the order of the distance 
between adjacent monomers a (the spatial resolution of 
the the harmonic chain). In comparison, temperature 
factors from x-ray crystallography correspond to RMS 
fluctuation of less than lA (a/4). A more realistic non- 
Gaussian polymer backbone would increase the resolu- 
tion of the model. Such a non-Gaussian model could still 
be treated using the variational formalism. Still retain- 
ing the Gaussian description of the chain, multi-body 
forces arising from neglected degrees of freedom would 
also tend to increase the rigidity of the native state and 
are a likely contributor. Recent off-lattice simulations 
have shown that effective multi-body forces sharpen in- 
terfaces between folded and unfolded regions and increase 
activation barriers E5 As-jdicussed in [I], contiguous se- 
quence approximationsEjilj also give inherently sharp in- 
terfaces and larger activation barriers even for pairwise 
additive forces since the conformations of the residues 
are restricted to be either folded or unfolded. (The con- 
tiguous sequence approximations with multi-body forces 
would give even higher barriers.) If multi-body forces are 
indeed important, the completely folded or unfolded con- 
formations assumed in these approximations may well be 
more appropriate than expected from simulations of lat- 
tice and off-lattice models with only pairwise forces. Still, 
the interfaces are probably not as cleanly sharp as envi- 
sioned in the contigous sequence approximations. Us- 
ing the calculated prefactor at the transition midpoint, 
fco ~ 8/iS -1 , and the measured folding rate of 100s _1 
gives a barrier of approximately llfcsTf. This value is 
greater than both the activation barrier from the varia- 
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tional theory (AF^ as 5ksTf) and the double sequence 
approximation calculated in [I] (AF^ w SfcsTf). Multi- 
body forces may sharpen the distinction between folded 
and unfolded regions in the transition state ensemble to 
give an activation barrier consistent with the prefactor 
calculation at the transition midpoint. 

Alternatively, it is possible the discrepancy at Tf arises 
because of neglected dynamical effects that could slow 
the chain dynamics, reducing the folding rate prefactor. 
In this paper, the internal friction is interpreted as a cor- 
rection to account for the activated isomerization reac- 
tion of the backbone. As mentioned in the introduction, 
there are other possible contributors to this friction. For 
example, non-native interactions are presently neglected 
in the model. These interactions contribute to the mem- 
ory friction and therefore the internal viscosity. They can 
be treated through mode coupling calculations similar to 
those we have done for a random heteropolymera (see 
also Refs. |36|-jB8l). A more careful treatment of the sol- 
vent may also be necessary for an accurate description 
of the chain dynamics. Detailed calculations of the reac- 
tion dynamics for isomerization transition in di-alanine 
(a much less complicated system) suggest the solvent may 
be an essential component of even the local reaction coor- 
dinate.!^ Furthermore, contact formation dynamics may 
also be slowed because contact formation excludes the 
solvent, thereby adding an additional activation barrier 
and timescale into the dynamics. Supporting this possi- 
bility, we note that simulations of minimal protein fold- 
ing models with desolvation barriers in the pair potential 
bring simulations into qualitative agreement with pres- 
sure studies on folding rates o 

Since we expect some of these neglected internal fric- 
tion contributions to influence protein folding dynamics, 
we have been careful in this paper to use "bare" parame- 
ters (and investigate other reasonable estimates) instead 
of treating the them as fitting parameters in the the- 
ory. This is an attempt to avoid missing a signal that 
such new dynamical considerations are indeed important 
by using inappropriate effective parameters. The folding 
rate alone, however, is not enough to determine the accu- 
racy of the calculated prefactor, since the barrier height 
may not yet be accurately enough detetermined. 

Consequently, these proposed effects are perhaps best 
studied at this point using more direct experimental 
probes such as intrachain fluorescence quenching mea- 
surements. Fluorescence quenching is a .powerful tool to 
study the average structure of ensernblesElo and has re- 
cently been used as a probe for experiments on single, 
molecules that allow distributions to be addressed. i3c3 
Motivated primarily by the relevance to protein folding 
dynamics, some fluorescence quenching experiments have 
focused jqii less complicated systems such as short pep- 
tides oo These measurements of simple polymer pair 
dynamics can guide the development of improved the- 
oretical models helpful to understanding protein fold- 
ing kinetics. This is particularly evident in the present 
model, since the expression for the effective mobility 



[Eq.(^||)] is formally related to appropriations to the in- 
trachain quenching rate in polymers ExE3 

Experiments monitoring the intrachain contact for- 
mation through fluorescence quenching suggest that the 
dynamical issues mentioned above play a role in the 
pair contact dynamics as well. Fitting data to effec- 
tive models for the pair dynamics, the apparent dif- 
fusion coefficient comes out to be 10-20 times smaller 
than one would-expect Ao-cn the diffusion coefficients of 
the free probes.ooEj'OLj In those models, the poly- 
mer is treated as a single pair connected by an entropic 
spring controlling the mean square fluctuation of chain 
enda£3 where the effective diffusion coefficient of the pair, 
D c s = 2Dq, is a short time approximation to the full 
dynamics including all the monomers. E3 Haas has sug- 
gested that the internal friction of the chain may explain 
the descrepancy between D e s and the fitted value,0 but 
no detailed calculation, has been done. This is a problem 
we wish to return to.E3 Doi's arguments suggest that the 
short time properties may dominated We note that the 
inclusion of internal viscosity from the bond flipping rate 
as presented here (7//7opHl0) reduces the effective short 
time diffusion coefficiently to a value D c g w -Do/2, nar- 
rowing the discrepancy to within a factor of about 2-5. 
Since the validity of this short time approximation can 
be questioned (e.g., the rate depends on the quenching 
volumeBES) , we feel a solution resolution to this puzzle 
requires a more complete calculation. 

The primary experimental technique used to study 
the structure of the transition state-,ensemble is the <f>- 
value analysis developed by Fersht£22l By comparing the 
change in the folding rate to the change in stability upon 
mutation, one can infer the extent to which that site is 
structured in the transition state ensemble. To interpret 
these experiments, one assumes that the mutation al- 
ters the kinetics through a change in the free energy pro- 
file, leaving the prefactor unchanged. Sometimes, how- 
ever, the kinetics and thermodynamics are seen not to 
be directly coupled. This may well be a sign of frustra- 
tion from non-stabilizing native contacts or non-native 
contacts influencing the kinetics of the transition.liil It 
is also possible however that these mutations affect the 
chain dynamics by modifying the internal friction at that 
residue. It is easy to extend the present calculation to 
indicate how the internal friction could be exploited to 
find "dynamical mutants" in which a residue is replaced 
by one with higher dihedral barriers (i.e., a residue with 
larger local internal friction) perturbing the prefactor, 
without altering the stability. Proline isomerization cou- 
pled to folding are an example of this, but other non- 
natural amino acids may be still more promising probes 



of local reaction coordinates 



Mutations of the n-src 



loop of the SH3 domain have very small effect on the 
equilibrium constant (preventing the <f>- value analysis), 
but still impact the folding-xate and thus may be exam- 
ples of dynamical mutants.E£3 More extreme alterations 
to the sequence have been necessary such as lengthening 
the loop length in order to characterize this region of the 
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transition state ensembleJ123 A better understanding of 



the chain dynamics should help suggest additional and 
more direct probes of similar situations using dynamical 
mutants. 



ACKNOWLEDGMENTS 

This work has been supported by NIH Grant No. PHS 
2 R01 GM44557. 



APPENDIX A: 



In this appendix, we outline the calculation of the na- 



tive density correlations (pi(t)pj(0))o given by Eq.(p6[). 

Since Hq is harmonic, the Green's function for the 
diffusion operator is known and has the form of a 
multi|0£ariable Gaussian containing all the pair correla- 
tions E3 The native density correlations only depend 
on two monomer positions for a fixed i and j, so that 
Eq.(p7j) can be simplified by introducing delta functions 
Jdvidr2 6(ti — Yi)8(r2 — Tj) and integrating over all 
monomer variables {r}. This reduces Eq.(p7|) to a double 
integral over ri and r2, and the reduced Green's function 
contains only the monomer correlations of the monomer 
pair, as can be expected for a Gaussian process. 

It is perhaps easier to make use this Gaussian property 
from the beginning, and obtain the same result as using 
the full Green's function as described above. 

Consider a fixed i and j pair. First we define some 
notation. We define the following 2-component vectors: 



R 



rj -(f) 



s = 



(Al) 



where, t and t' denote two times, and e.g., s, = (n)o and 
rf is the average and native position vector of the i th 
monomer, respectively. Let x(^~ be the symmetric 

2x2 matrix of monomer correlations x^^') 1 = (8R(t)- 
SK(t') T )o/a 2 , with 8R(t) = R(t) - S: 



x(\t-t' 



Gij{\t-t'\) 



Gij(\t-t'\) 



(A2) 



To construct x(\t ~ f\): we have denoted the equal time 
correlation function by the static correlations Gu and 
Gjj, and used symmetry of the correlations Gij(t) — 
Gji(t). 

Since this is a Gaussian process, the correlation func- 
tion (pi{t)pj(t))o can be written as 

(Pi(t)pj(t')} = [ pt(T i )p j {T j )P(ST i ,6T j ;\t--e\), (A3) 

where 



P(*r i ,&- J ;|t-t'|) 



2na [ 



■(det X (|i-t'|)) 3 / 2 



x exp 



_3 

2a 



2 8R T ■ X (\t - f\) ■ STL 



(A4) 



In the rest of the derivation, we suppress the time de- 
pendence to simplify notation. Substituting Eq.(^) for 
p{vi)p{Yj) and shifting the integration variables to SR 
(R = <5R + S) leads to 



(pi(t)pj(t))o = ^2(det X ) 3/2 exp 



3c/ 
' 2a 2 



(R N - S) 2 



x / exp 

' Sri,Srj 



~5R-jS -cTR+^J SR 
2cr a z 



where 



x' = X + al 



I = 



1 
1 



J = R y 



(A5) 



(A6) 



(A7) 



Performing this integral by completing the square in the 
exponent gives 



(Pi(t)Pj(t))o 



det\ 



1 3/2 



det x' 
x exp 



exp 



3c/ 



3(a N ) 2 
2a 2 



2a 2 



(A8) 



To simplify this equation, it is convenient to re-write 



X as 



x' = xM 



(A9) 



where 



M={I + a N X^) 

l + a N G 2l a N G J3 (|i-i'|) 
a N G. ti (|i - l + « N G y - 

Then we can express X as 

X^ = M-\-\ 

Now since 

I = M~ X M 

= M- 1 (I + a N X - 1 ) 
= M _1 +a N M" 1 x" 1 
= M- 1 +a N X '- 1 , 



we can express x as 



X 



(A10) 



(All) 



(A12) 
(A13) 
(A14) 
(A15) 



(A16) 



Consequent ly, t he combined argument the exponential 
factors in Eq.rtAq) is 
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3(a N ) 2 
2a 2 



1 



(J J • X • J - -W J 



3a N 
" 2a 2 



J J • M • J. 



(A17) 

The determinate factor can be simplified using 

clct x' = det x det M. (A18) 
Thus, Eq.(|A8|) becomes 

3a N 



(Pi(t)pj(t))o = (dctAf)- 3 / 2 exp 



which is Eq. 



2a 2 



(A19) 
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TABLE I. 
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ture, Tf. The parameters are: I = 5a 
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TABLE II. Calculated kinetic parameters at two sta- 
bilities In A = ~P(Fm — Fg) where Ajv and Fg are the 
free energy of the native and globule states: contact energy 
eo/fcsT, barrier height AF' /ksT, corresponding folding pref- 
actor ko = \k*\/2tt, and folding time r = Other parame- 



ters are fixed: I = 5 a (g — O.i 
71/70 = 9.1. 



), a = 0.5, Do = 3 x 10- 6 cm 2 /s, 



A? 2 ° x 10 6 [cm 2 /s] 




71/70 


ko [^s _1 ] 


3 


3 


9.1 


7.6 


3 


10 


2.5 


8.7 


9.1 


3 


28 


18 


9.1 


10 


8.2 


23 



TABLE III. Internal friction and prefactor (from TS2) 
for A-repressor for different values of the monomer diffu- 
sion coefficient in water D^ 2 ° and isomerization rate r^*. 
D^ 2 ° i=n9.1 x 10~ 6 cm 2 /s is the diffusion coefficient of free 
alanineH and r^h_= 10 ns -1 is the calculated isomerization 
rate for di-alanineo . The parameters in the first line are con- 
sidered the best estimates. The persistence length is I = 5 a 
(g — 0.8), the width of the Gaussian measure of the native 
density is a = 0.5, and the temperature is the folding transi- 
tion temperature. 
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FIG. 1. Mode amplitude plotted vs. sequence index for 
the four modes with the slowest relaxation rates; from the 
slowest mode: solid, dashed, long dashed, dotted. (7/ = 0.) 
The constraints correspond to stationary points on the free 
energy surface for A-repressor with persistence length I — 5a 
(g = 0.8): (a) unconstrained, (b) Globule, (c) TSi, (d) TS 3 , 
(e) TS4,, (f) Native. 
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FIG. 2. Relaxation rates A p vs. mode index for different 
constraints corresponding to stationary points on the free en- 
ergy surface for A-repressor with persistence length I — 5a 
(g = 0.8): unconstrained (C = 0), TSi, TSa, and Native (N) 
are indicated in the plot. 
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of internal viscosity: 7//7O = and 7//7O = 10 (indicated 
on the plot). The constraints correspond to A-repressor with 
persistence length I — 5a (g — 0.8): unconstrained (dashed 
lines) and TS2 (solid lines). 
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FIG. 5. Local reaction coordinates for the transition states 
TS1-TS4 of the A-repressor average folding route at Tf. Nor- 
malization of the reaction coordinate is set to ^2 — 1. 
Parameters are the same as in Fig.^. 
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FIG. 4. Native density pi evaluated at the native and 
globule minima (dashed) and the transition states TS1-TS4 
(solid) of the average folding route for A-repressor. The per- 
sistence length of the chain is I = 5 a (g = 0.8) and the width 
defining p; is set by a — 0.5. 
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FIG. 6. Prefactor ko = \k*\/2ty in units of a = 3Do/a 2 as 
a function of 71 /70 for the transition states TS1-TS4 of the 
A-repressor average folding route at Tf. Parameters are the 
same as in Fig.0. 
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FIG. 7. 
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Prefactor fco is plotted as a function of 



the inverse monomer diffusion coefficient D n 



The internal friction is set equal to the fixed value 
Ji/ksTf = 9.1/D^ 2 °. The vertical dotted line corresponds 
-Do = A? 2 ° = 3 x 1CT 6 cm 2 /s (and 77/70 = 9.1). Parameters 



are the same as in Fig.^. 



5.0, 



4.0 



3.0 



CO 

o 



2.0 



o 



1.0 



0.0 [ 




0.2 



TS1 



TS4 



0.4 



0.6 



0.8 



a 



N 



FIG. 8. Prefactor ko — \k*\/2tt vs. the Gaussian width of 
the native density a N . The persistence length of the chain is 
I — 5a (g — 0.8) and the internal friction is 71/70 = 0.91. 
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